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PACS 05. 40. -a - Fluctuation phenomena, random processes, noise, and Brownian motion 
PACS 47.27.eb - Turbulent flows; statistical theories and models 

Abstract. - We reveal a phase transition with decreasing viscosity f at 1/ — Vc > in one- 
dimensional decaying Burgers turbulence with a power-law correlated random profile of Gaussian- 
distributed initial velocities v{x,0)v{x' ,0) ~ |a: — a;'|"^. The low-viscosity phase exhibits non- 
Gaussian one-point probabihty density of velocities, continuously dependent on v, reflecting a 
spontaneous one step replica symmetry breaking (RSB) in the associated statistical mechanics 
problem. We obtain the low orders cumulants analytically. Our results, which are checked nu- 
merically, are based on combining insights in the mechanism of the freezing transition in random 
logarithmic potentials with an extension of duality relations discovered recently in Random Ma- 
trix Theory. They are essentially non mean-field in nature as also demonstrated by the shock size 
distribution computed numerically and different from the short range correlated Kida model. We 
also provide some insights for the finite viscosity behaviour of velocities in the latter model. 



Analysis of the solutions of the Burgers equation 9tv + 
(vV)v ~ z^V^v -|- /(x, i) with either random forcing 
/(x, t) 7^ or unforced but subject to random initial con- 
ditions v(x, t = 0) 7^ (the latter case being frequently 
referred to as the "decaying Burgers turbulence" (dBt)) 
attracted a considerable interest in the last two decades, 
see [1] for an accessible introduction and a detailed lit- 
erature overview. The problem appears as an important 
reference model not only in fluid dynamics, but also in 
such diverse physical contexts as statistics of growing in- 
terfaces [2] , statistical mechanics of systems with quenched 
disorder [3-7], and formation of large scale structures in 
cosmology [8]. In particular, the cosmological applica- 
tions stimulated interest in dBt for vanishing viscosity 
I' — > and scale-free power-law random initial conditions 
v(x, 0)v(x', 0) ~ |x — x'l^"^^ at large distance, see e.g [9] 
and references in [1]. The latter are more properly defined 
via the Fourier-transformed value, the mean initial energy 
spectrum Eo(k) - fc''-iv(k, 0)v(-k, 0) - |k|", d being 
the space dimension. The model can be mapped, via the 



Cole-Hopf transformation [1], to the statistical mechan- 
ics of a particle in a d dimensional random potential at 
temperature T = 2v. Although intensive analytical and 
numerical work resulted in a good qualitative understand- 
ing of general features of emerging velocity structures, de- 
tailed analytical results are mainly available in d = 1 for 
the following particular cases: (i) the white-noise initial 
velocity [6, 10] (n = 0) (ii) the Kida case of short range 
(SR) potentials {n = 2) [11] (iii) Brownian motion initial 
velocity {n = —2) [12]. It is believed that the general 
case n < 1 leads to a self-similar evolution with persis- 
tence of the fc" spectrum i?t>o(fc) at small k (called per- 
sistence of large eddies in turbulence [13], or equivalently 
long-range (LR) fixed points in the RG study of the par- 
ticle model [6]), while n > 1 leads to Kida's SR spec- 
trum, with however, a crossover region 2 > n > 1 where 
both types of behaviours can be found depending on the 
scale [13,14]. Interestingly, it seems that the limiting case 
n ~ 1 has not been studied in any detail although as we 
show below it is in many respects special. The goal of our 
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Letter is to fill in that gap. 

To achieve this we develop a method allowing one to 
get analytical insights into the velocity statistics of dBt 
for the "marginal" case n = 1. We concentrate on analyz- 
ing the simplest quantity, namely the probability density 
hmction (p.d.f.) 7^(1') of a single point velocity v = v{x, t), 
at any viscosity i> > 0. Our main finding is that the shape 
of the function V{v) experiences an abrupt change at a 
finite viscosity value v ^ Vc > 0- Namely, for v > Vc 
the probability density is Gaussian, which is the natural 
and expected result, but 7^(u) ceases to retain the Gaus- 
sian shape in the low-viscosity phase, i.e. everywhere for 
Q < V < Vc- We reveal that such a change is in fact a 
manifestation of the so-called freezing phase transition in 
the associated statistical mechanics problem of a single 
particle in random logarithmically correlated landscape. 
The transition of such a type was discovered long ago in 
mean-field type models [15], see [16] for recent activity, 
and is identified with the simplest (one-step) pattern of 
spontaneous replica symmetry breaking (RSB). More re- 
cently, accumulated evidences pointed towards existence 
of a similar transition in finite dimensions [17-19], al- 
though well-known challenges of extending RSB pattern 
beyond the mean-field level remained outstanding. In 
the present study we show how to incorporate the one- 
step RSB into the finite-dimensional calculation by fur- 
ther developing and adopting to the situation a method 
suggested in [5]. Such a procedure when combined with 
recent progress in Random Matrix Theory (RMT) yields 
highly nontrivial predictions for the velocity cumulants 
in the low viscosity phase. Comparing them to numeri- 
cal results on the particle model which solves the Burgers 
equation shows a good agreement and provides a rather 
convincing support to validity of the method. The fact 
that our results are essentially of non mean-field nature 
is also demonstrated by the shock size distribution, which 
we compute numerically and find different from the short 
range correlated Kida model, itself weh described by a 
mean field one step RSB ansatz. Finally, we get some 
analytical insights into the finite viscosity behaviour of 
velocities in the delta-correlated case by extending [5] . 

Our starting point is the standard mapping [1] of the 
dBt problem to the equilibrium statistical mechanics of a 
single classical particle in a potential V{x) at the effective 
temperature T played by the viscosity: T = 1v. The 
initial velocity profile v{x^ 0) in this approach is related 
to V{x) as u(a;,0) = ^y(a;). To ensure n ~ \ power- 
law spatial decay of the initial conditions we choose the 
random potential Yix) to be Gaussian with the two-point 
function given by 

V{x)V{x') = -21n[|x-a;'|/-^], e<\x~x'\<L (1) 

where i ^ 1 and e <C 1 are the infrared and ultravi- 
olet cutoff scales, correspondingly. We further assume 
V{x)V(x') = 2\TLL/e for \x-x'\ < e and V(x)Vix') = 
for [a; — a:' I > L. 



The solution to the unforced Burgers equation for t > 
with chosen initial conditions is given by v(x,t) = dxV, 
where V{x,t) = —T\nZv{x,t) is the "renormalized po- 
tential" (effective free energy functional) corresponding to 

the Hamiltonian 'H{y;x) = -I- F(j/), with the parti- 

tion function defined as 

Zy{x,t) = ef^ j ^ _^cxp--H(2/;x). (2) 

where /3 = and the e-factor is chosen to facilitate the 
comparison with our earlier work [19]. 

To this end we note that in the language of statis- 
tical mechanics the velocity p.d.f. is given by V{v) = 
(5 (u + J -< y ^t) where we introduced the thermal aver- 
age ^ O >T= Zy^J -^O(y)exp-H(y;0)/r and set 
x = in view of the translational invariance of the dis- 
order. To understand better thermodynamics of our sys- 
tem and the nature of the anticipated freezing transition 
it turns out to be instructive to consider also a different 
object: Vyiy) = -< 5{Y — y) interpreted as the aver- 
aged p.d.f. of the coordinate of a particle equilibrated 
at a given temperature T in the random energy land- 
scape H{y;0). At r ^ the thermal average is obvi- 
ously dominated by the deepest minimum of the land- 
scape whose position j/mm fluctuates from one realization 
of disorder to the other. This mechanism immediately im- 
plies for velocity p.d.f. in zero viscosity limit the relation 

riv)\T=o^trYivt)\T=o- 

The disorder averaging procedure for VyiY) can be per- 
formed via the standard replica trick after representing 
Zy^ = Zy~^\n^Q and using the Gaussian nature of the 
random potential V{y). Employing Eq.(IT]) yields the rela- 
tion 

where j ~ (3'^ > and we have defined for 1 < ?? < I/7 

(4) 

with 5„(-7) = n'-" [r(l-J7)/r(l-7)] being the fa- 
mous Selberg integral. For finite integer n > 1 and 
—7 = A > the above expression is nothing else but 
the mean density of the so-called A-Hcrmitc ensemble of 
RMT introduced by Dumitriu and Edelman (DE) [20]. 
Although a closed-form expression for the eigenvalue den- 
sity for that ensemble does not seem to be available yet, 
DE developed analytic tools to compute a few lower mo- 
ments of that density for any integer n > 0. We no- 
ticed that their result can be recast in terms of cumu- 
lants as < zf >^ 1 — 7(n — 1), and for integer q > 1, 

< z^ >fj _..y= — 7(71— l)(7n — l)P2g(?^,7) in terms of poly- 
nomials of n and 7, which read to low orders Pi{n, 7) = 1, 
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PQ{n,-i) = -2 + 7(5n-2). P8(«,7) = 6 + 267 + 67^ - 
47771 — 477^71 + 567^7j^. They satisfy a remarkable du- 
ality relation 'y'^^'^P2q{nj, ^) = P2q{n,'^), and since their 
order does not grow too fast, we could obtain them by 
direct evaluation of < z^"^ >n.-7 for (small) integer 77 and 
7 < 0, up to order 2q ~ 16 (we also used DE's algorithm 
up to 2q = 12 as a check). Performing the analytical 
continuation 7i — > and < 7 < 1 (where the correspond- 
ing RMT-like integrals are still convergent) and assuming 
So (—7) = 1 we obtained the lower nonvanishing moments 
M2q = J VY{Y)Y'^''dY up to 2q = 16, i.e. in term of their 
cumulants Ca, =-< 7/^« )-f,= -(Ti)?7P2g(0, 7), yielding: 
C2=t{T + T-i) , Ci = -e, Ce = 2^3 (r + T-i) 

Cg^-t^ [26 + 6(r2 + r-2)] 

Clo = [300 (T + T-^) + 24 (T^ + T'^)] 

and similar but longer expressions for C2q, q = 6, 7, 8. The 
second cumulant agrees with the exact relation ^ _ 
t^v"^ + Tt valid at any T due to statistical translational in- 
variance of the random potential V{x). The main feature 
apparent from the above (and proved in full generality) 
is that all the cumulants (and hence the whole function 
Vyiyy) are invariant with respect to the duality transfor- 
mation T — > 1/T. Duality of this kind was first discovered 
in our work on logarithmically-correlated random poten- 
tuals [19] where it was conjectured that it implies freez- 
ing transition at the self-dual point Tc = 1: the self-dual 
functions retain down to zero temperature the shape they 
acquired at the critical point T = Tc- Here we thus pre- 
dict that the whole probability distribution VyiY) freezes 
at T = 1 providing a vivid picture of what freezing entails. 
If this scenario were correct, the values of the above cu- 
mulants evaluated at T = 1 should immediately provide, 
in view of the discussed zero-temperature correspondence, 
the cumulants of the velocity p.d.f. in zero viscosity limit: 
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Fig. 1: Numerical evaluation of and in the inviscid limit 
r = for M = 2^°, 2", 2^* compared to the prediction (O 
at t = I (averaged over 10*^ samples). Small oscillations are 
observed in when Ma ~ 1 and the periodic boundary con- 
ditions cannot be neglected, and disappear when Ma ^ 1. 



A 

?^ 
Y 




2 -TC 



1 



72 (5) 



and more generally w^9|i/=o = t ^'^C2q\T=i- Due to the ^ ^ j^^g 



above exact relation it also predicts 



Fig. 2: The top solid line is the analytical prediction (from 
C2 in the text) for -< ^t, and the bottom solid line for 
t^Ti^ = ^ y y'^ for t — 1 from (|14|) . Circles are simulations with 
M = 2^*, triangles are simulations with a — 1/1024, 



(2 - T) in M = 2^* (averaged over 5 x 10* samples) 



the whole low-T phase, as recovered below. Note that at 
T = all positive integer moments of v^/t are integersl 

The above cumulants are the nontrivial predictions of 
the developed approach, and are now tested against a nu- 
merical study of the Hamiltonian H(j/;a;). We discretize 
y = ai and x = aj with i,j = 1, ..M and study H{i;j) = 

° — I" where Vi are gaussian log-correlated vari- 
ables. Generally, a sample of M correlated variables has 
high numerical cost (^ or operations), a major 
simplification occurs when the correlation matrix is cir- 
culant and direct diagonalization via Fast Fourier Trans- 
form is possible (M log M operations) Our results are ob- 
tained using the log-circular ensemble [18, 19] of gaussian 
variables with correlations ViVi' ~ —2 ln2| sin( ^'''^' ^ )|, 



i 7^ i', and ~ 2\nM. For this ensemble it is pos- 
sible to generate a sample up to M ^ 2^° points. This 
ensemble, for L = Ma 3> 1 reproduces the model on 
the line ([1]) for jx — ^ L. We compute in the dis- 
crete model ^ (i — j) which in the limit a — > at 
fixed L = Ma ^ 1 identifies to the velocity of the con- 
tinuum model v{x = aj,t = 1) with a cutoff e ^ a. This 
yields the result at all times for the continuum model, us- 
ing the exact rescaling v{x,t)\e^L 
which shows self-similarity in the limit e <^ t <^ L 

As a first check of the freezing scenario, one can see 
in Fig. [1] that the prediction ([5]) is well obeyed if one 



7i^(7^'l)le/Vt,L/Vt 
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takes into account finite size corrections as a — >■ 0. It is 
apparent from this result that the velocity distribution is 
non-gaussian in the inviscid limit T = 0. A second check, 
see Fig. [2l is the tem perature dependence of both -< >t 
and v"^ = < y >~'ip . In the high-T phase the agreement 
with the analytical prediction is perfect. For T < 1 the 
convergence becomes much slower, but in agreement with 
the freezing predictions, e.g. that -< should freeze 

at a constant value for T < 1. 

After satisfying ourselves with the validity of the freez- 
ing scenario, we arc going to briefly demonstrate how such 
a transition can be recovered from the replica calcula- 
tion, which eventually yields the expression for the veloc- 
ity cumulants for any viscosity below critical. To illustrate 
the method we consider the simplest object, the partition 
function moments 
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V2TrTt 



1 



(6) 

where n is the replica number eventually set to zero. In 
doing this we adopt to the continuum model the scheme 
of incorporating one-step RSB mechanisms proposed in [5] 
for the simplest case of discrete Random Energy Model 
without spatial correlations. The basic idea behind this 
scheme is that for T < Tc and < n < 1 the configurations 
which give the leading-order contributions to the above 
integral are obtained by grouping n replica indices into 
k = n/m groups of m replica each, and assuming that 
all coordinates yi for the replica indices zi, . . . , inside 
the same group are "frozen" around the common value, 
i.e. approximately equal: y^j^ ~ « . . . ~ y^^ . More 
precisely, they are allowed to fluctuate within a distance 
of order of the small-scale cutoff e around their common 
centre of mass j/i, . At the same time k coordinates 

of the centres of masses of different groups play the role of 
new effective degrees of freedom and can take any values. 
Integrating out the "frozen" coordinates yields the factor 
of the order of e"^'^, so that Eq.® is replaced with 



7" c"'' +n~k/^ 
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where the combinatorial factor Cn,m ~ nl/klml'^ takes 
into account the number of ways we can built the groups. 
At this stage we can perform the disorder average in the 
standard way using Eq.(IT]) and find that the above ex- 
pression is proportional to the to— dependent large factor 
exp — [n -f- ^) Ine] . The parameter m is then found 
from extremizing (in fact, minimizing) this factor, which 
selects TO = T as long as T < Tc = 1. Performing 
the calculation to the end one reproduces, up to a con- 
stant shift in free energy, the expressions of the moments 
Zy, < n < 1 above and below Tc = 1, precisely those 



following from the freezing scenario for logarithmic mod- 
els [18,19]. The details of this and other calculations out- 
lined in the paper will be given elsewhere [22], here we 
give only very brief account. Following the method of [19] 
we first find exactly arbitrary complex moment of the nor- 
malized partition function Zy — Zv(27r)~2~r(l — /3^). In 
the high-temperature phase /3 < 1 they are given by 



2^2 



(8) 



where s is arbitrary complex and we use the same conven- 
tion for the generalized Barnes function Gp{x) as discussed 
in [19]. Note that the factor = L" ^ (which amounts to a 
convolution with a gaussian of large variance 2/3^ In L ^ 1) 
ensures the convexity of the moments and makes the prob- 
lem of restoring the corresponding probability distribu- 
tion of Zv well-defined. This expression, combined with 
the conjecture of freezing of all self-dual functions allows 
to restore after due manipulations the partition function 
moments below the freezing temperature. The results co- 
incide with the direct replica calculation outlined above. 

To find the shape of the velocity p.d.f. ^^{v) for T > 
it is convenient to exploit the generating function G{q) = 
In [1 — iqdxV{x^t)] which can be further calculated via a 
variant of the replica trick G{q) = lim„^o ^(W^n('z) — 
W^n(O)) where for integer n > 



k=0 



(9) 



Following the same steps as before we arrive at the identity 



W„(<z) 



\i=i 



n,— 7 

(10) 



To continue to ti = we exploit the relation 




Aw 



27r 



[r + i^/XwY (11) 



where r is an arbitrary parameter. Although it was origi- 
nally proved in the RMT context assuming integer n> 1 
and A > (see Theorem 4.1 in [20] ) we conjecture it 
can be analytically continued beyond the original domain 
as long as all integrals make sense. Applying Eg. pT]) for 
A = — 7 with < 7 < 1 and performing the replica limit 
n — > in Ea.(|10p we find after due manipulations the ve- 
locity probability density Vlv). It turns out to be a simple 
Gaussian with zero mean and variance v'^ = 1/Tt. 

The idea of freezing suggests that the above result 
should be valid as long as T > 1. To calculate the gener- 
ating function G{q) for the velocity moments in the low- 
temperature phase r < Tc = 1 we again have to employ 
the same 1-step RSB scheme as described earlier in the 
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Fig. 3: Non-Gaussian character of the rescaled velocity v = 
v/Vv^ below the freezing temperature, as shown by the differ- 
ence between their cumulative distributions. Simulations are 
performed over 10^ samples of M = 2^* points and a = 0.002. 
Circles are data at T = 0, Squares are data at T = 0.5 < Tc. 



paper. Up to factors tending to unity in the replica limit 
n — > we arrive at the relation 



Wniq)^ {Y[(l + tqVt 



tzij ) , m = T 

f k—^,—'ym^ — — l 

(12) 

Finally, we notice that one can perform the replica limit 
ri — >■ by exploiting a powerful duality relation for 
A— Hcrmite RMT ensemble discovered in [21] for k, m pos- 
itive integer, A > and any complex s: 



— mA 



-feA 



(13) 



\i=i 



k.X 



\l=l 



m. A 



We conjecture that the relation remains valid if continued 
to A = — 1 and furthermore to < k,m < 1. This al- 
lows to perform straightforwardly the n — )■ limit leading 
to the expression of the velocity p.d.f Viv) = lim,„_j.T < 
S{v + ^) >„i,-i as eigenvalue density in the DE ensem- 
ble. At T = it indeed identifies with ([3]) which confirms 
the freezing scenario. At any < T < 1 the velocity cu- 
mulants are thus again given by the above polynomials 
= lim„,^T < >m,-i= -(1-T)2p2g(r, 1), which 
we computed up to 2q = 16. Hence we find 



V\t<1 



-(2-r), 



|T<1 = 



1 



[l^Tf (14) 



which fully agrees with the zero-viscosity limit and 
matches high-temperature phase moments at the transi- 
tion point Tc = 1- These results are in agreement with 
numerics in Fig. [U and show that the velocity p.d.f. V{v) 
is non-Gaussian everywhere in the low-viscosity phase, as 



seen in Fig. [31 The shape is consistent with a negative 
kurtosis and the difference increases at low temperature. 

We now put our results in context, compare with the SR 
Kida behaviour (obtained here by replacing ViVj = crSij), 
and discuss shocks. Let us recall the expected scaling 
for the LR fixed points n < 1, either from FRG argu- 
ments for the particle problem [6], or from the Burgers 
literature [1, 13] (although there n = 1 was not specifi- 
cally discussed): self similarity holds with v{x,t) =iniaw 
ti^^v{xt~'^/'^) and energy decay E(t) ~ l/t~^+'', with 
C = 4/(3 + n) the Flory value. For n < 1 the energy 
exponent for the particle, 9 ~ 2{( — 1) > 0, with the 
effect that the size of the viscous layer of the shocks van- 
ishes when rescaled by the average shock distance. The 
case n ~ 1 studied here clearly belongs to this LR family 
with = 0, C = 1 and indeed velocity and energy decay 
respectively as v ~ l/Vi and E{t) ^ 1/t [23]. This is 
distinct from the Kida SR behaviour E{t) - l/t{\iit)^^^ . 
We expect that the relative shock width/f^/^ remains con- 
stant and T dependent, consistent with 9 — marginal- 
ity (Fig. 4, inset). The (one point) shock size distri- 
bution at T = has been computed numerically, with 
S = v{x~^,t) — v{x~ ,t). The SR case reproduces correctly 
Kida's exact resuh [11] P{S) - S'^^'e"^"^' with t = -1. 
However, the present model yields a different distribution, 
as seen in FigH) Examination of the cumulative distri- 
bution (not shown here) reveals that the exponent r is 
consistent with t « which would be the prediction of 
the general conjecture [24] r = 2 — -gq:^- With D = 
for decaying Burgers it predicts r = for the LR fixed 
points n < 1, which also recovers the exact result r = 1/2 
for n = [6, 10]. It shows that despite its one step RSB 
character the present model is quite different from mean 
field behaviour (which in the one step case always yields 
P{S) - Se-^'^' [25]). 

Finally, one can apply our method to the SR Kida prob- 
lem, hereby extending to T > the results of [11] and 
adding to some of the observations in [5] . Considering the 

with (T 



aS, 



above discrete model with ViVj 
and t = l/a^, one finds for K t/p < AP 



,2 - 



-(( 

rMn(27rt/p) 



1/2 



T) 



2ttT 



0(1), 



e^(15) 



where p = m/T is solution of the saddle point equation 
ap^ = l+ln(27ri/p). The freezing corresponds to the range 
< m < 1 and the condition m ~ 1 yields the above 
value of tc for a given temperature T. For t <^ tc one 
recovers the T = decay law of Kida, while for t > tc the 
energy has fully decayed and v = 0. Eq thus describes 
what happens in-between, as the system effectively heats 
up: for t < tc the system is in the glass phase (i.e. here 
Tc ^ (Ini)^^/^ obtained setting m = 1), tc being the time 
at which the thermal width of the shocks ^ Tt has reached 
the distance between shocks ~ '^(ei)^^^- This behaviour 
is shown in Fig. ([S]). Computing for that model Wn{q) 
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Fig. 4: Distribution of the rescaled shock sizes s = S/S at T — 
0. Continuous line, Kida's prediction p(s) — ^se . Open 
circles: SR Kida discrete model with 1/a = 256. Filled circles, 
squares: logarithmic model, for 1/a = 64, 1024 respectively 
with M = 2^^. Inset: plot of a typical velocity field v{x,t = 1) 
for the logarithmic model with M = 2^"^ and 1/a = 1024 and 
2u = T = 0, 1/3, 10. The jumps at T = are the shocks. 



and setting v ~ v/ ^Jtp we obtain 



ln(l — iqv) 



In I 



dy 



(1 + iqyY 



SO that = {\ — m) and v"^ = — 2to(1 — to) and 



where S± 



|~|-2-2mg-CV2 



27rr(l - to) Si + S± +2cos{7rm)S+S. 



/c 



oo 




27rr(l-?ri) 



"2 (^±1) . The asymptotics is 
3""^/^ at large v and P{0) is finite. 



Extension of the present theory to velocity- velocity cor- 
relations, some moments of shock size distribution, and 
comparison with the infinite-dimensional limit will be 
given elsewhere [22] . The transition unveiled here in d = 1 
does extend to any dimension d for n = 1 (i.e. in- 
verse square spatial decay initial velocity correlation) with 
Tc = 2vc = \fd [17]. We expect that qualitatively the be- 
haviour of velocity will be similar to the d = \ case. 
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